In this study, we have steroid information of 31 male patients at three different time points after trauma (<1h, 4-12h, 24-48h). Only 26 have complete steroid information for the three time points. We had also immunological information for these patients but it was taken out from the study. Clinical information such as age, MODS, LOS, ISS, Drugs, Fluids.., is also included. 34 healthy controls are also included for reference.
Two analysis with different steroid’s concentration have been performed given we had two different LOQ values given. The first measurements taken with old LOQ and the new revalidated LOQ.
This is the one with the REVALIDATED LOQ
#Loading packages
#install_github("kassambara/easyGgplot2")
require(readxl)
require(reshape2)
require(plyr)
require(lubridate)
require(stringr)
require(ggrepel)
require(data.table)
require(ggpubr)
require(readr)
require(tidyverse)
require(dplyr)
require(knitr)
require(htmltools)
require(factoextra)
require(Metrics)
require(e1071)
require(ggplot2)
require(gridExtra)
require(brglm)
require(ggstatsplot)
require(compareGroups)
require(knitr)
require(kableExtra)
require(gtsummary)
require(dunn.test)
require(broom)
require(reshape2)
require("patchwork")
require(gam)
require(mgcv)
require(ComplexHeatmap)
#this might need github too
#packageurl <- "https://cran.r-project.org/src/contrib/Archive/repr/repr_1.0.2.tar.gz"
#install.packages(packageurl, repos=NULL, type="source")
set.seed(132)
File with LOQ concentrations (lowest concentration at which the analyte can not only be reliably detected).It includes the parameter we used at the start, and at which the original analysis was developed (original) and the corrected/revalidated one, where many datapoints were filtered as a result.
Again this file contains all analysis performed with REVALIDATED LOQ
Steroids <- read_xlsx("SteroidsLow.xlsx")
Steroids <- Steroids[complete.cases(Steroids),c(2,3,4)]
names(Steroids) <- c("variable", "LOQ_Original", "LOQ_Revalidation")
34 healthy controls with steroid information. Standard ranges were also recovered and measurements outside those ranges were taken out.
#outliers per steroid.
#saveRDS(Names, "names.rds")
RealHealthy <- function(HealthyLong){
Names <- data.frame(variable = unique(filter(HealthyLong, Descriptor == "Steroid")$variable))
Names["LimitLow"] <- c("30","10","1","0.7","2","0.1","2","0.6","1","0.09","0.3","2","0","0.07","0","1", rep("0",7)) #said 1 in progesterone (second to last but then all readings would be off)
Names["LimitHigh"] <- c("100","650","2.6","3.9","13","0.7","40","4","7","0.6","40","35","10","2.4","10","20",rep("100",7)) #much less in170HP but then off too * also 11KA4 very high numbers in general - muh than trauma and reference
Names
a <- full_join(HealthyLong,Names)
a$LimitLow <- as.numeric(a$LimitLow)
a$LimitHigh <- as.numeric(a$LimitHigh)
a1 <- a %>%
mutate(TakeOut = value >= (LimitLow) & value <= (LimitHigh)) #take out ten percent of readings? -- Yes
HealthyAll <- a1 %>%
filter( TakeOut == "TRUE") %>%
select(-c("LimitLow","LimitHigh", "TakeOut"))
return(HealthyAll)
}
Healthy <- read_delim("FinalHealthyPat.csv",
";", escape_double = FALSE, locale = locale(decimal_mark = ",",
grouping_mark = "."), trim_ws = TRUE) %>%
select(-Age) #Take out Age for healthy - should be age matched either way.
names(Healthy)<- c("ID", "Label", "Cortisone", "Cortisol", "11-KTestosterone","11-KA4", "11β-OHA4","11β-OHT", "Corticosterone", "11-deoxycortisol", "Androstenedione", "11-deoxycorticosterone", "Testosterone", "DHEA", "17-hydroxyprogesterone", "Dihydrotestosterone", "Progesterone", "DHEAS", "Cortisol/Cortisone", "DHEA/DHEAS", "Cortisol/DHEAS", "Cortisol/DHEA", "DHEA/Testosterone", "Testosterone/Cortisol", "Aldosterone") #distinguish between trauma and healthy control
HealthyLong <- reshape2::melt(Healthy, id.vars=c("ID","Label"))
HealthyLong <- HealthyLong %>%
mutate(Time = "HC",
Descriptor = "Steroid")
HealthyLong <- RealHealthy(HealthyLong) #LOQ feeds out of this one so all fine
HealthyLongLOQ <- full_join(HealthyLong,Steroids) %>%
replace_na(list(LOQ_Original = 0, LOQ_Revalidation = 0)) %>%
filter(value >= LOQ_Revalidation)
HealthyWideLOQ <- HealthyLongLOQ %>%
select(ID, Label, variable, value, Descriptor) %>%
pivot_wider(names_from = variable, values_from = value)
#take out healthy patients that do not follow normal ranges (??)
31 trauma patients with cytokine, immunological and clincal information at three different time points.
(Lesson: pivot longer will put problems with types and classes of columns but melt wont).
Trauma <- read_csv("FinalDatasetAprilwCytokines.csv") #non healthy patients
names(Trauma)[1] <- "ID"
TraumaLong <- reshape2::melt(Trauma, id.vars=c(1))
TakeOut <- c(grep("^[*]",TraumaLong$value),which(TraumaLong$value %in% c("ND","OOR <"))) #Take out spurious values
TraumaLong <- TraumaLong[-c(TakeOut),] %>%
drop_na() %>%
mutate(Label = "Trauma")
TraumaLong[,c("variable","Time") ] <- str_split_fixed(TraumaLong$variable, "-", 2)
TraumaLong$Time <- gsub(" ","",TraumaLong$Time)
TraumaLong$Time <- as.factor(TraumaLong$Time)
levels(TraumaLong$Time) <- c("Clinical","T1","T2","T3")
'%ni%' <- Negate('%in%')
TraumaLong <- TraumaLong %>%
filter(variable %ni% c("IL1ra","IL6", "IL8","IL10","GCSF", "MCP1", "TNFa","IL12","IL17", "Lymphocyte (106/L)","Monocyte (106/L)","Neutro:Lymph ratio","Neutrophil (106/L)","WBC (106/L)","IG (106/L)")) %>%
mutate(Descriptor = ifelse(Time == "Clinical", "Clinical", "Steroid"))
TraumaLong$variable <- as.factor(TraumaLong$variable)
levels(TraumaLong$variable) <-c("11-deoxycortisol","11-KA4","11-KTestosterone","11β-OHA4","11β-OHT","17-hydroxyprogesterone", "Androstenedione","Age","Alcohol", "Aldosterone","Anaesthetic","Corticosterone","Cortisol",
"Cortisol/Cortisone","Cortisol/DHEA","Cortisol/DHEAS","Cortisone","DHEA","DHEA/DHEAS","DHEA/Testosterone","DHEAS",
"Dihydrotestosterone","11-deoxycorticosterone","Fluids","GCS","GoldenHour","ISS","MODS", "Lactate","LOS","Mechanisms","NISS","Outcome","Drug","Progesterone","Sex","Testosterone","Testosterone/Cortisol",
"T1","T2","T3")
TraumaWide <- TraumaLong %>%
select(- c(Label, Descriptor)) %>%
pivot_wider(values_from = value, names_from= c(Time,variable))
TraumaLongLOQ <- full_join(TraumaLong,Steroids) %>%
replace_na(list(LOQ_Original = 0, LOQ_Revalidation = 0)) %>%
filter(value >= LOQ_Revalidation)
names(TraumaLongLOQ)[1] <- "ID"
TraumaWideLOQ <- TraumaLongLOQ %>%
select(ID, Label, Time, variable, value) %>%
pivot_wider(names_from = c(Time,variable), values_from = value)
CreateTablesTP <- function(TraumaWide,TraumaWideLOQ, TimePointOut1,TimePointOut2, TimePointIn ) {
TraumaT1 <- TraumaWide %>%
select(-c(contains(TimePointOut1), contains(TimePointOut2), "ID")) %>%
#drop_na() %>%
mutate_at(vars(starts_with(TimePointIn), "Clinical_Age", "Clinical_GCS", "Clinical_ISS", "Clinical_NISS", "Clinical_LOS", "Clinical_Lactate", "Clinical_GoldenHour" ),funs(as.numeric))
names(TraumaT1) <- str_remove_all(names(TraumaT1), TimePointIn)
TableTraumaT1 <- TraumaT1 %>%
tbl_summary(
type = list(all_numeric() ~ "continuous",
all_logical() ~ "categorical") ,
missing_text = "(Missing)"
) %>%
bold_labels() %>%
modify_header(label = "**Variable**")
#
TraumaT1LOQ <- TraumaWideLOQ %>%
select(-c(contains(TimePointOut1), contains(TimePointOut2), "ID")) %>%
#drop_na() %>%
mutate_at(vars(starts_with(TimePointIn), "Clinical_Age", "Clinical_GCS", "Clinical_ISS", "Clinical_NISS", "Clinical_LOS", "Clinical_Lactate", "Clinical_GoldenHour" ),funs(as.numeric))
names(TraumaT1LOQ) <- str_remove_all(names(TraumaT1LOQ), TimePointIn)
TableTraumaT1LOQ <- TraumaT1LOQ %>%
tbl_summary(
type = list(all_numeric() ~ "continuous",
all_logical() ~ "categorical") ,
missing_text = "(Missing)"
) %>%
bold_labels() %>%
modify_header(label = "**Variable**")
return(list(TableTraumaT1,TableTraumaT1LOQ ))
}
TP1 <- CreateTablesTP(TraumaWide,TraumaWideLOQ, "T2_","T3_", "T1_" )
TP2 <- CreateTablesTP(TraumaWide,TraumaWideLOQ, "T1_","T3_", "T2_" )
TP3 <- CreateTablesTP(TraumaWide,TraumaWideLOQ, "T2_","T1_", "T3_" )
##########################
HealthyTable <- Healthy %>%
#drop_na() %>%
select(-c(ID, Label)) %>%
tbl_summary(
type = list(all_numeric() ~ "continuous",
all_logical() ~ "categorical") ,
missing_text = "(Missing)"
) %>%
bold_labels() %>%
modify_header(label = "**Variable**")
HealthyTableLOQ <- HealthyWideLOQ %>%
#drop_na() %>%
select(-c(ID, Label, Descriptor)) %>%
tbl_summary(
type = list(all_numeric() ~ "continuous",
all_logical() ~ "categorical") ,
missing_text = "(Missing)"
) %>%
bold_labels() %>%
modify_header(label = "**Variable**")
tbl_merge(tbls = list(HealthyTableLOQ, TP1[[2]], TP2[[2]], TP3[[2]]), tab_spanner = c("**Healthy Controls**", "**Time Point <1h**", "**Time Point 4-12h**", "**Time Point 24-48h**") )
| Variable | Healthy Controls | Time Point <1h | Time Point 4-12h | Time Point 24-48h |
|---|---|---|---|---|
| N = 341 | N = 312 | N = 312 | N = 312 | |
| Cortisone | 62 (53, 70) | 53 (42, 60) | 57 (50, 70) | 34 (22, 45) |
| (Missing) | 3 | 2 | 5 | |
| Cortisol | 273 (210, 335) | 458 (386, 573) | 510 (325, 605) | 214 (114, 321) |
| (Missing) | 2 | 5 | ||
| 11-KTestosterone | 1.50 (1.16, 1.60) | 1.4 (0.8, 2.4) | 10 (1, 12) | 0.8 (0.7, 11.8) |
| (Missing) | 16 | 2 | 5 | |
| 11-KA4 | 3.20 (2.65, 3.45) | 1.8 (1.4, 2.8) | 1.69 (1.12, 3.01) | 1.50 (1.10, 3.97) |
| (Missing) | 23 | 2 | 5 | |
| 11β-OHA4 | 8.45 (6.32, 10.07) | 16 (12, 22) | 12 (8, 21) | 9.6 (5.2, 12.4) |
| (Missing) | 12 | 2 | 8 | |
| 11β-OHT | 0.600 (0.600, 0.625) | 0.92 (0.78, 1.23) | 0.79 (0.78, 1.61) | 0.6886 (0.6886, 0.6886) |
| (Missing) | 27 | 9 | 25 | 30 |
| Corticosterone | 7 (5, 15) | 61 (49, 88) | 39 (15, 76) | 7 (3, 15) |
| (Missing) | 2 | 8 | ||
| 11-deoxycortisol | 1.10 (0.80, 1.68) | 4.5 (3.0, 7.6) | 5.5 (3.0, 7.7) | 1.17 (0.66, 2.58) |
| (Missing) | 8 | 1 | 6 | 15 |
| Androstenedione | 2.97 (2.62, 4.32) | 3.98 (3.24, 5.69) | 3.20 (2.15, 5.45) | 3.12 (2.09, 4.13) |
| (Missing) | 3 | 12 | ||
| 11-deoxycorticosterone | 0.30 (0.30, 0.35) | 0.74 (0.40, 1.64) | 0.95 (0.81, 2.32) | 0.77 (0.56, 0.98) |
| (Missing) | 27 | 9 | 13 | 29 |
| Testosterone | 14.5 (11.5, 17.4) | 11.9 (8.7, 13.9) | 3.83 (2.63, 5.75) | 1.7 (0.9, 4.5) |
| (Missing) | 2 | 7 | ||
| DHEA | 19 (13, 21) | 32 (24, 54) | 18 (14, 33) | 8 (3, 15) |
| (Missing) | 5 | 2 | 5 | |
| 17-hydroxyprogesterone | 3.25 (1.83, 4.07) | 5 (4, 7) | 4 (2, 11) | 1.37 (0.93, 2.69) |
| (Missing) | 3 | 15 | ||
| Dihydrotestosterone | 1.60 (1.15, 1.87) | 1.08 (0.65, 1.29) | 0.60 (0.60, 0.77) | 0.60 (0.60, 0.60) |
| (Missing) | 8 | 2 | 6 | |
| Progesterone | 0.4000 (0.4000, 0.4000) | 0.76 (0.60, 1.26) | 1.04 (0.80, 1.93) | 0.344 (0.339, 0.346) |
| (Missing) | 31 | 3 | 12 | 28 |
| DHEAS | 9.3 (6.4, 11.1) | 13 (8, 20) | 8 (6, 18) | 14 (10, 18) |
| (Missing) | 2 | 5 | ||
| Cortisol/Cortisone | 4.53 (3.88, 6.16) | 8.87 (8.03, 10.14) | 9.3 (5.8, 11.5) | 5.58 (4.57, 7.35) |
| (Missing) | 2 | 5 | ||
| DHEA/DHEAS | 2.24 (1.56, 3.74) | 3.00 (1.73, 3.87) | 2.04 (1.33, 3.54) | 0.77 (0.22, 1.55) |
| (Missing) | 2 | 5 | ||
| Cortisol/DHEAS | 30 (22, 44) | 37 (26, 55) | 55 (26, 86) | 15 (8, 25) |
| (Missing) | 1 | 2 | 5 | |
| Cortisol/DHEA | 13.6 (10.4, 17.1) | 13 (10, 19) | 110 (56, 230) | 22 (13, 43) |
| (Missing) | 2 | 5 | ||
| DHEA/Testosterone | 1.34 (0.92, 2.02) | 3.19 (2.03, 4.91) | 4.2 (2.8, 10.6) | 3 (2, 9) |
| (Missing) | 2 | 5 | ||
| Testosterone/Cortisol | 0.049 (0.039, 0.065) | 0.024 (0.019, 0.033) | 0.009 (0.004, 0.018) | 0.010 (0.004, 0.027) |
| (Missing) | 2 | 5 | ||
| Aldosterone | 0.60 (0.60, 0.60) | 0.74 (0.60, 1.07) | 0.70 (0.60, 1.05) | 0.60 (0.60, 0.60) |
| (Missing) | 2 | 5 | ||
| Label | ||||
| Trauma | 31 (100%) | 31 (100%) | 31 (100%) | |
| Clinical_Age | 26 (22, 30) | 26 (22, 30) | 26 (22, 30) | |
| Clinical_Sex | ||||
| Male | 31 (100%) | 31 (100%) | 31 (100%) | |
| Clinical_GCS | 15.00 (14.00, 15.00) | 15.00 (14.00, 15.00) | 15.00 (14.00, 15.00) | |
| Clinical_ISS | 16 (10, 21) | 16 (10, 21) | 16 (10, 21) | |
| Clinical_NISS | 22 (16, 29) | 22 (16, 29) | 22 (16, 29) | |
| Clinical_Fluids | 9 (29%) | 9 (29%) | 9 (29%) | |
| Clinical_Anaesthetic | 12 (39%) | 12 (39%) | 12 (39%) | |
| Clinical_Drug | ||||
| Cocaine | 1 (3.2%) | 1 (3.2%) | 1 (3.2%) | |
| Fentynl | 1 (3.2%) | 1 (3.2%) | 1 (3.2%) | |
| Ketamine | 6 (19%) | 6 (19%) | 6 (19%) | |
| Morphine | 10 (32%) | 10 (32%) | 10 (32%) | |
| NO | 11 (35%) | 11 (35%) | 11 (35%) | |
| Paracetemol | 1 (3.2%) | 1 (3.2%) | 1 (3.2%) | |
| Rocuromium | 1 (3.2%) | 1 (3.2%) | 1 (3.2%) | |
| Clinical_MODS | 13 (42%) | 13 (42%) | 13 (42%) | |
| Clinical_Mechanisms | ||||
| Assault | 1 (3.2%) | 1 (3.2%) | 1 (3.2%) | |
| Assault and Penetrating | 4 (13%) | 4 (13%) | 4 (13%) | |
| Blunt Impact | 2 (6.5%) | 2 (6.5%) | 2 (6.5%) | |
| Fall From Height | 1 (3.2%) | 1 (3.2%) | 1 (3.2%) | |
| Penetrating | 8 (26%) | 8 (26%) | 8 (26%) | |
| RTC | 15 (48%) | 15 (48%) | 15 (48%) | |
| Clinical_Alcohol | 5 (16%) | 5 (16%) | 5 (16%) | |
| Clinical_Outcome | ||||
| Alive | 31 (100%) | 31 (100%) | 31 (100%) | |
| Clinical_LOS | 10 (6, 20) | 10 (6, 20) | 10 (6, 20) | |
| Clinical_Lactate | 3.20 (2.60, 4.78) | 3.20 (2.60, 4.78) | 3.20 (2.60, 4.78) | |
| Clinical_GoldenHour | 36 (26, 46) | 36 (26, 46) | 36 (26, 46) | |
| Clinical_T1 | ||||
| Afternoon | 13 (42%) | 13 (42%) | 13 (42%) | |
| Evening | 7 (23%) | 7 (23%) | 7 (23%) | |
| Morning | 4 (13%) | 4 (13%) | 4 (13%) | |
| Night | 7 (23%) | 7 (23%) | 7 (23%) | |
| Clinical_T2 | ||||
| Afternoon | 12 (41%) | 12 (41%) | 12 (41%) | |
| Evening | 4 (14%) | 4 (14%) | 4 (14%) | |
| Morning | 12 (41%) | 12 (41%) | 12 (41%) | |
| Night | 1 (3.4%) | 1 (3.4%) | 1 (3.4%) | |
| (Missing) | 2 | 2 | 2 | |
| Clinical_T3 | ||||
| Afternoon | 2 (7.7%) | 2 (7.7%) | 2 (7.7%) | |
| Morning | 24 (92%) | 24 (92%) | 24 (92%) | |
| (Missing) | 5 | 5 | 5 | |
|
1
Median (IQR)
2
n (%); Median (IQR)
|
||||
We have more missing values in revalidated than original, because this new revalidation has caused many measurements to become invalid (missing). To start with we had 5 missing patient values in T3 and 2 for T2.
AllDataLOQLong <- rbind(TraumaLongLOQ, HealthyLongLOQ)
AllDataLong <- rbind(TraumaLong, HealthyLong)
General steroid concentration changes at the three different time points
(k3[["Progesterone"]] | k3[["17-hydroxyprogesterone"]]) /
(k3[["11-deoxycortisol"]]| k3[["Cortisol"]]) /
(k3[["Cortisone"]] | k3[["Cortisol/Cortisone"]])
Figure.1
k3[["11-deoxycorticosterone"]] | k3[["Corticosterone"]] | k3[["Aldosterone"]]
Figure.2
(k3[["DHEA"]] | k3[["DHEAS"]] | k3[["DHEA/DHEAS"]]) /
(k3[["Androstenedione"]] | k3[["Testosterone"]] | k3[["Dihydrotestosterone"]])
Figure.3
(k3[["11β-OHA4"]] | k3[["11β-OHT"]]) /
(k3[["11-KA4"]] | k3[["11-KTestosterone"]] )
Figure.4
Study of steroid concentration at the first time point (<1h). Use of Generalized Additive Models (GAMs) to extract trajectories through time.
wrap_plots(pQuanta) + plot_layout(guides = "collect")
Golden Hour GAM
Original plots with patients above/below normal healthy concentrations. GAMs tries to build two different models (around blue dots and orange dots). As you can see, those steroids without model are those where the separation between orange and blue is practically non existent (e.g DHEA/DHEAS)
wrap_plots(pQuant2a) + plot_layout(guides = "collect")
Golden Hour All
Mineralocorticoids <- c("11-deoxycorticosterone","Corticosterone","Aldosterone")
GlucocorticoidPrec <- c("Progesterone","17-hydroxyprogesterone","11-deoxycortisol")
Glucocorticoids <- c("Cortisol", "Cortisone" )
Androgens <- c("DHEA","DHEAS","Androstenedione")
OxygenatedAndro <- c("Testosterone", "Dihydrotestosterone","11-KA4","11-KTestosterone","11B-OHA4","11B-OHT")
Ratios <- c("Cortisol/Cortisone", "Cortisol/DHEA","Cortisol/DHEAS","DHEA/DHEAS","DHEA/Testosterone","Testosterone/Cortisol")
aa <- list(Mineralocorticoids,GlucocorticoidPrec, Glucocorticoids,Androgens, OxygenatedAndro, Ratios)
names(aa) <- c("Mineralocorticoids","GlucocorticoidPrec", "Glucocorticoids","Androgens", "OxygenatedAndro", "Ratios")
ll <- data.frame(aa) %>%
pivot_longer(everything()) %>%
unique() %>%
as.data.frame()
https://github.com/kassambara/ggpubr/issues/65 p value adjust
# it was working before and now not - issues with names.
scale_this <- function(x){
(x - mean(x, na.rm=TRUE)) / sd(x, na.rm=TRUE)
}
Heatmapping <- function(ChosenData,ll){
AverageHC <- ChosenData %>%
filter(Time =="HC") %>%
group_by(variable) %>%
mutate(MeanValue = mean(as.numeric(value))) %>%
select(variable, MeanValue) %>%
unique()
#AverageHC$variable <- str_replace(AverageHC$variable, "11β-OHA4", "11B-OHA4")
#AverageHC$variable <- str_replace(AverageHC$variable, "11β-OHT", "11B-OHT")
MixingHM <- inner_join(ChosenData, AverageHC) #adds a column with mean HC value
FinalDataHM <- MixingHM %>%
filter( Descriptor %in% c("Steroid"), Time !="HC")
FinalDataHM$variable <- str_replace(FinalDataHM$variable, "11β-OHA4", "11B-OHA4")
FinalDataHM$variable <- str_replace(FinalDataHM$variable, "11β-OHT", "11B-OHT")
FinalDataHM$variable <- droplevels(as.factor(FinalDataHM$variable))
FinalDataHM$value <- as.numeric(FinalDataHM$value)
FinalDataHM <- FinalDataHM %>%
mutate(valueWHC = (value-MeanValue)) %>%
group_by(variable) %>%
mutate(value_scaled = scale_this(valueWHC)) %>%
select(ID, variable, value_scaled, Time) %>%
pivot_wider(id_cols = ID, names_from = c(variable, Time), values_from = value_scaled)
Clinical <- ChosenData %>%
filter( Descriptor != c("Steroid"), Time !="HC") %>%
select(ID, variable, value, Time) %>%
pivot_wider(id_cols = ID, names_from = c(variable, Time), values_from = value) %>%
as.data.frame()
FinalDataHMID <- FinalDataHM$ID
FinalDataHMMatrix <- as.matrix(FinalDataHM %>%
select(-c(ID, '11B-OHT_T3','11B-OHT_T2'))) ## too many missing values
rownames(FinalDataHMMatrix) = FinalDataHMID
ll3 <- data.frame(str_split(colnames(FinalDataHMMatrix), "_", simplify = TRUE))
names(ll3) <- c("value", "time")
ll2 <- inner_join(ll3, ll)
annotation_col = data.frame(
Time = ll2$time,
Type = ll2$name
#Mean = ll2a$Mean
)
rownames(annotation_col) = colnames(FinalDataHMMatrix)
annotation_row = data.frame(
Age= as.numeric(Clinical$Age_Clinical),
LOS = as.numeric(Clinical$LOS_Clinical),
ISS = as.numeric(Clinical$ISS_Clinical),
Fluids = as.factor(Clinical$Fluids_Clinical),
Anaesthesia = as.factor(Clinical$Anaesthetic_Clinical)
)
rownames(annotation_row) = Clinical$ID
ann_colors = list(
Type = c(Androgens = "light blue", GlucocorticoidPrec = "yellow", Glucocorticoids = "orange", Mineralocorticoids = "springgreen4",OxygenatedAndro= "dark blue", Ratios ="snow2"),
Time = c(T1 = "green", T2 = "orange", T3 = "red"),
Anaesthesia = c(NO = "lawngreen", YES = "red2"),
Fluids = c(NO = "lawngreen", YES = "red2")#,
#Age = c("white", "firebrick"),
#ISS = c("white", "firebrick"),
#LOS = c("white", "firebrick")
)
#cutree_cols = 5,cutree_rows = 3,
Heatmap <- pheatmap(FinalDataHMMatrix, annotation_col = annotation_col, annotation_row = annotation_row, cutree_cols = 5,cutree_rows = 3,annotation_colors = ann_colors, fontsize_col = 9, labels_row =paste0("Patient.", 1:31))#, fontsize_row = 0
# annotation_colors = ann_colors)
return(list(FinalDataHM, ll,Clinical, Heatmap ))
}
Once again, main difference with original LOQ data is the amount of missing values (grey cells) [ Had to delete ‘11B-OHT_T3’,‘11B-OHT_T2’ as had only one value left]
AldNo <- AllDataLOQLong %>%
filter(variable != "Aldosterone")
Heatmapping(AldNo,ll)[[4]]
Heatmap with NO LOQ change
pdf("HeatmapAll.pdf", 17, 11)
Heatmapping(AldNo,ll)[[4]] #too many missing values!
dev.off()
## quartz_off_screen
## 2
# it was working before and now not - issues with names.
scale_this <- function(x){
(x - mean(x, na.rm=TRUE)) / sd(x, na.rm=TRUE)
}
Average <- AldNo %>%
filter(Time != "Clinical") %>%
group_by(variable) %>%
mutate(value_scaled = scale_this(as.numeric(value))) %>%
ungroup() %>%
group_by(variable,Time) %>%
mutate(MeanValue = mean(as.numeric(value_scaled))) %>%
select(variable, MeanValue) %>%
unique()
Average$variable <- str_replace(Average$variable, "11β-OHA4", "11B-OHA4")
Average$variable <- str_replace(Average$variable, "11β-OHT", "11B-OHT")
AverageMat <- pivot_wider(Average,names_from = Time, values_from =MeanValue )
FinalDataHMID <- AverageMat$variable
AverageMatFinal <- AverageMat %>% as.data.frame() %>%
dplyr::select(c("HC","T1", "T2", "T3"))
FinalDataHMMatrix <- as.matrix(AverageMatFinal)
rownames(FinalDataHMMatrix) <- FinalDataHMID
ll2 <- ll[match(rownames(FinalDataHMMatrix), ll$value),] %>%
filter(value != "Aldosterone")
annotation_row = data.frame(
Type = ll2$name
)
rownames(annotation_row) = rownames(FinalDataHMMatrix)
ann_colors = list(
Type = c(Androgens = "light blue", GlucocorticoidPrec = "yellow", Glucocorticoids = "orange", Mineralocorticoids = "springgreen4",OxygenatedAndro= "dark blue", Ratios ="snow2")
)
######
AveragePlot <- pheatmap(FinalDataHMMatrix, annotation_row = annotation_row, annotation_colors = ann_colors, fontsize_col = 9,cluster_cols=FALSE)#,
AveragePlot
Heatmap average
pdf("HeatmapAverageConce.pdf", 7, 5)
print(AveragePlot)
dev.off()
## quartz_off_screen
## 2